a=zeros(50,50);
for i=1:50
    a(i,i)=12;
end
for i=1:49
    a(i+1,i)=-2;
    a(i,i+1)=-2;
end
for i=1:48
    a(i+2,i)=-1;
    a(i,i+2)=-1;
end
b=ones(50,1);
jac=inv(diag(diag(a)))*(-1*(tril(a,-1)+triu(a,1)));
gs=inv(tril(a))*(-1*triu(a,1));
pbj_jac=max(abs(eig(jac)))
pbj_gs=max(abs(eig(gs)))
x1=[];
x2=[];
x=zeros(50,1);

eps = 0.00001

for k=1:50
    x=jac*x+inv(diag(diag(a)))*b;  
    x1=[x1,x];
    if k > 1
        dx = abs(x1(1:50, k) - x1(1:50, k - 1));
        if dx < eps
            k
            break;
        end
    end
end

subplot(1,2,1)
plot(x1)
hold on
x=zeros(50,1);
for k=1:5
    x=gs*x+inv(tril(a))*b;  
    x2=[x2,x];
    if k > 1
        dx = abs(x2(1:50, k) - x2(1:50, k - 1));
        if dx < eps
            k
            break;
        end
    end
end

subplot(1,2,2)
plot(x2)